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Summary of the report. We present a fast direct solver for the volume scattering 
problem of the Helmholtz equation. The algorithm is faster than existing methods. 
Moreover, discretization for our method is much simpler and more accurate than that 
for finite difference, finite elements, and integral equations. 

Jacques Hadamard's work on ill-posedness put us in a box of solving well-posed prob- 
lems of preferably small condition number. In reformulating elliptic problems, such as 
the scattering problem for the Helmholtz or Maxwell, by integral equations (IEs) the 
price we have to pay is the complexity in solution representations and their discretiza- 
tions, and cost of computation (20 or more, instead of 2 to 4, points per wavelength) 
for a fast direct solver (FDS). Well-posedness lead us, for example, to reformulating the 
first order Maxwell equations as second order elliptic PDEs, with symmetries breached 
and balances of physical quantities disturbed. Additional difficulties will have to be met 
when the second order equations are formulated by IEs, and as the latter are discretized 
and solved numerically. 

There is another all-embracing box next to that of Hadamard's. The more capable 
we become of solving problems, the more we seek challenging and interesting problems, 
such as a scattering problem. In the present report we propose a different approach to 
solution of elliptic problems in a compactly supported domain D with inhomogeneous 
medium, referred to here as the volume scattering problem (VSP). Our method never 
solves a well-posed or ill-posed problem. Instead it first solves a very simple but not 
posed problem: PDEs without boundary conditions. 

We will construct the total wave solution space (TWSS) for the given PDEs, or more 
precisely, we will construct the null space of the homogenous PDEs (with zero RHS and of 
variable coefficients), subject to no boundary or any other conditions. Thus the TWSS 
consists of general solutions of the PDEs. Constructing the null space seems a non- 
scattering problem but it may be the easiest way to account for global communictions of 
the linear system, or global multiple scatterings of a scattering problem. It is only after 
the TWSS is constructed in D that the data, or the incident wave, is incorporated to 
obtain a specific solution - the specific total wave corresponding to the specified incident 
wave. 

Our method will give rise to a fast direct solver for the elliptic problem in D. It 
does not support an iterative solver. The TWSS will be constructed recursively on a 
(quadtree for example, in 2-D) hierarchy of domain decomposition, with the TWSS first 
constructed in each bottom level subdomains. Merging the total waves in the subdomains 
to those of their parents will end up with the TWSS constructed efficiently for the entire 
scatterer D. 

The method will be presented in the scattering context and language, but the prin- 
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ciples extend to a general boundary value problem for elliptic PDEs, in particular the 
Maxwell equations. There are 8 first order equations for the 6 unknown E and B for lin- 
ear material. We can construct the null space of the linear equations in D first, and deal 
with the boundary conditions latter. Another major benefit of the total wave approach 
is that it greatly simplifies discretization. 

In this report, a total wave refers to a nontrivial solution of the homogenous PDEs 
with variable coefficients in domain D. 

1 Introduction 

The subject of this report is solution of scattering problem for the Helmholtz equation. 
There are two standard types of scattering problems. One is for the inhomogeneous 
medium inside a domain and is referred to here as the volume scattering problem (VSP). 
The other is for an impenetrable scatterer (such as a perfect conductor for Maxwell), 
or a penetrable scatterer with constant coefficients inside, and is referred to as surface 
scattering problem (SSP). When formulated as integral equations (IEs), VSP is related 
to the Lippmann-Schwinger IE, whereas SSP is related to boundary IEs. 

These problems can also be solved by finite element or boundary element methods. 

There are fast direct solvers (FDSs) and iterative solver. A FDS is desirable when 
the condition number of the problem is not small, which occurs, for example, when the 
problem is near resonance. They are also more efficient for a scattering calculation with 
multiple incident waves typically required for an inverse scattering problem. 

For IEs or finite element methods, discretization has never been made robust, or even 
easy. For a second kind IE, with all its underlying benefits in conditioning and reduction 
of dimensionalities, the discretization problem is even more evident. The use of layer 
potentials makes their discretization extremely unwieldy and difficult. For example, a 
quadrature to integrate the Coulomb potential 1/r near y, with r = \x — y\ and with y a 
fixed point on a patch of smooth surface, against smooth functions such as polynomials, 
is not easy to design due to the strong and non-trivial influence of the curvatures on 
1/r near the source y. Warping a smooth function we get a smooth function. Warping 
a singular function, we'd better prepared to reap the whirlwinds. So far we have not 
looked at the frightening situations when these singular kernels meet corners and edges 
on surfaces or inside inhomogeneous media. 

The desirable analytic properties and attractiveness of IE formulation for the scat- 
tering problems, or for any elliptic PDEs, or for Maxwell equations, become difficult to 
exploit the moment they meet discretization. This is because after discretization we have 
to deal with individual and standing alone poles and dipoles. A pole will not interact 
well with other nearby poles or singularities unless they are on the same patch of smooth 
surface. Unfortunately, for many interesting applications, different parts of the surface, 
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or multiple inclusions, may get very close. On the fly, designing a quadrature will be 
more difficult than resampling. 

We propose a solution method for the scattering problem, also for other elliptic PDEs 
and the Maxwell equations, which does not require discretization of IEs or PDEs. It 
requires what we call "sampling" of the original PDEs (as in a collocation method) and 
their solutions. We don't call it discretization in the sense that the collocation method 
may be regarded as sampling, as opposed to discretizing the differential equation. There 
is a distinction between sampling and discretization. Sampling requires not much of 
brain, whereas discretization requires too much of it. 

Our method does not reformulate the original PDEs as IEs. It works with the PDEs 
directly. It barely solves those PDEs. Certainly it never solves a well or ill-posed problem 
for the PDEs. It solves a not posed problem for the PDEs. It solves the PDEs without 
boundary or other conditions. For scattering problems or other elliptic PDEs formulated 
as a boundary value problem, the boundary values will be processed only after the general 
solution space for the PDEs without boundary conditions are constructed. Constructing 
the general solution space is easier than solving boundary value problems. 

The method is a fast direct solver; it cannot be related to and does not support any 
iterative solver. For a SSP, with the surface not very concave or convoluted, our FDS is 
faster than existing FDSs. For VSPs, our FDS offers the same asymptotic complexity, 
with a big reduction on the constant if the inhomogeneous medium occupies a convex 
domain D, such as a square or triangle. 

Organization of the report: £f2]is a short and informal description of the method. $3]is 
a full and more formal description of the method. §4] provides formulations of the volume 
scattering problem and Green's third identity used as a projector on the boundary of a 
scatterer. §5] contains background information on layer potential representation for the 
interior and exterior projectors. 



2 Our method - Informal description 

We will present the method in the context of VSP for the Helmholtz; see $H for more 
details on VSP. For a given precision e > and incident wave Uq, the algorithm finds 
the unique solution to the scattering problem in a compactly supported inhomogeneous 
medium (a variable index of refraction n(x)) inside domain D. If v denotes the scattered 
wave, then u = uq + v is the total wave, which satisfies the homogeneous Helmholtz with 
variable coefficients 

Au + k 2 n 2 (x)u = 0, xeD (1) 
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2.1 Informal description 

Our method consists of three parts. 

Part I. For the prescribed precision e > 0, a complete set of solutions {uj,j = 1 : N} 
to the homogeneous Helmholtz ([1]) are efficiently constructed in D and made available 
on the boundary dD. The boundary lies in the free space. See $3]for further details on 
the size N and how the N solutions are constructed efficiently. 

Part II. Each solution Uj thus constructed, being a total wave, can and will be 
split into two parts on dD, the incident and scattered waves uoj and Vj. This can be 
accomplished with the third Green's identity used as a projector on L 2 (dD). Since the 
set of total waves {uj} is complete, any total wave, namely a solution of the homogeneous 
Helmholtz (jTJ), can be represented by linear combination of {uj, j = 1 : N} to precision e. 
Likewise, any incident wave can be represented by linear combination of {uqj, j = 1 : N} 
to precision e. In particular, our prescribed incident wave u can be expressed in terms 
of {u 0j J = l:N} 

N 

n o( x ) = ^^ c j u oj( x )i x £ 9D (2) 
i=i 

Solving (J2J) for the coefficients Cj, we obtain the scattered wave 

N 

v i x ) = CjVj(x), x E dD (3) 
i=i 

corresponding to Uq. At this point, we have obtained the scattered wave on dD and 
consequently also outside D. For many applications, such as inverse scattering by re- 
peatedly solving forward problems, the scattered wave outside the medium is all that we 
want. 

Part III. Now suppose we also want the scattered wave v inside D. It will be obtained 
efficiently by a downward recursive procedure along a hierarchical structure, such as a 
quadtree for 2-D domain D, which was also used in Part I to efficiently construct the 
total waves {uj, j = 1 : N} in the first place. 

2.2 Discussions 

The hierarchical mergings for Part I, and splittings for Part III are universally employed 
in a typical FDS, although they may not always be presented in a familiar language or 
structure. It will referred to as domain decomposition. 

Definition 2.1 Throughout the report, D is regarded as singly connected. A subdomain 
is always the result of partitioning D artificially for the domain decomposition. 
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The method can also be adopted for SSP and scattering problems in layered media, 
so that no layered Green's function is required which is necessary for the IE formulation. 
See §3.31 for more details on the extensions to SSP. 

The method solves no scattering problem in order to construct the TWSS. It never 
solve a well or ill-posed problem for PDEs. In fact one can largely roam in the null 
space of linear operators without encounter ill-posedness; see §|2] for further details. The 
modern concept of first or second kind IEs seems oxymoron to the primitive kernel 
hunter-gatherer, and null spaces of IEs or geometric resonances are his trophies to hang 
up on the walls, not to become his stumbling blocks. 

A total wave in a subdomain, such as a square, subject to no condition on the 
boundary, is only aware of the medium inside the domain. It has no knowledge of what 
is outside, in particular whether the medium is discontinuous over the boundary, and 
thus it is unaware of the corners of the domain, unless explicitly informed; see (fl9l) . In 
contrast, a scattered wave is a solution of an inhomogeneous Helmholtz, subject to out- 
going radiation conditions in the free space, in which the subdomain must be embedded 
to set up the scattering problem for the subdomain, on whose solution a scattered wave 
based FDS relies. Thus the scattered wave sees the manmade discontinuity across the 
boundary. It is aware of the corners and requires more points there to be represented. 

3 Detailed description of the Algorithm 

We first present in §3.1l the basic components required by our fast direct solver. We then 
describe our TWSS based FDS in §3.21 The algorithm is similar to those of [l]-[5] in 
data structure and complexity; Our method differs from theirs in what solution space to 
construct and how to merge solution spaces of subdomains; see Remark 13.41 

3.1 The Basic Components and Parameters 

For simplicity, we assume that the scatterer q is a smooth function, which vanishes 
smoothly outside a square domain D. In implementation, this requirement can be relaxed 
to include piecewise smooth scatterers with jumps in a fairly arbitrary, bounded domain 
D in two dimensions. 

A typical fast direct solver for the Helmholtz equation with large wave number k relies 
on domain decomposition of some sort [1]- [8]. For the Lippmann-Schwinger equation 
( 135]) . the square domain D is partitioned hierarchically into the balanced quadtree; again 
for simplicity we will not discuss adaptive partitioning until §3.31 

The size of a scattering problem is measured by the number of wavelengths in each 
linear dimension, and the number of points required to discretize D is therefore propor- 
tional to k 2 . Let N = 0(k 2 ) be the number of unknowns in the resulting linear system 
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of equations to be solved. It is well known, [T], [I], that a fast direct solver requires 
0(iV L5 ) = 0(k 3 ) flops to construct general solution space and additional 0(N log N) 
flops to obtain the specific scattered wave for the prescribed incident wave. Our ap- 
proach has the same complexity, with a considerable reduction on the constant. 

Domain decomposition can be carried out, in principle, before or after the discretiza- 
tion. Our approach will be able to cleanly separate the two issues in a straightforward 
way, and partition D before discretization. In contrast, existing implementations go 
the other way around, and have to deal with "subdomains" of mesh points of the dis- 
cretization. Linear algebra tricks combined with untidy local approximation steps are 
employed to copy with these "subdomains" and the communications among them. The 
artificial cuts and corners must still manifest themselves in the "subdomains of points". 
For the Lippmann-Schwinger equation, the points are monopoles or dipoles. Strategies 
have been designed and strifes directed to these issues, to try to mend or heal the cuts 
and wounds. 

Observation 3.1 For a prescribed precision and in a subdomain, typically a square or a 
rectangle obtained by merging two squares, the number of distinct solutions, whether the 
incident waves, scattered waves, or total waves, in the subdomain is proportional to the 
arclength of the subdomain, as measured by the number of wavelengths. This is because 
these solutions can be determined uniquely by their Dirichlet, or Neumann, or D&N data 
on the boundary of the subdomain. 

Definition 3.2 For a prescribed precision the finite number of distinct solutions form 
the solution space for the subdomain in question. 

Remark 3.3 For simplicity, we will say that the dimension of solution space for a square 
subdomain, of edge length L and boundary arclength AL, is 4L, instead of proportional 
to AL. Therefore, the dimension of solution space for a rectangle by merging two squares 
is 6L. 

A necessary step of a fast direct solver is to construct the solution space for the entire D 
efficiently in a hierarchical order. Once solution space is formed for each subdomain on 
the same level of the quadtree, two neighbor square subdomains on the level are merged 
together to form the solution space for the union of two squares, and so on. 

Remark 3.4 The main difference between existing methods and our approach is that (i) 
We merger the total waves whereas they merge the scattered waves (ii) They organize the 
scattered waves by the incident waves whereas we don't organize the total waves by the 
incident waves or any other waves (Hi) We merge the total waves by simple continuity 
conditions whereas they merge the scattered waves via the multiple scattering process 
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among subdomains. (iv) We deal directly with solutions whereas they deal with operator 
such as the scattering matrix or the Dirichlet to Neumann map. 

Consequently, we never solve a scattering problem in merging; they solve an interesting 
and difficult scattering problem for each sub domain and in every step of merging. 

Definition 3.5 A total wave in D is a non-trivial solution of the homogeneous, variable 
coefficient Helmholtz equation l[3l\) in D, subject to no boundary conditions. A total 
wave in a subdomain is a non-trivial solution of (3l\) in that subdomain, subject to no 
boundary conditions, other than natural extension to the outside by continuity conditions 
on the function and its normal derivative on the boundary. 

As such, it has a desirable, simple property: A total wave in a subdomain is utterly 
oblivious to the underlying multiple scattering process among the subdomains. The 
total waves are much easier to construct and merge to build up the entire solution space 
for D. 

Observation 3.6 As two neighbor square subdomains merge, the two solution spaces, 
each of dimension AL, should merge to the solution space for the rectangle and of dimen- 
sion 6L . Indeed, the 2 continuity conditions (on function and its normal derivative) on 
the common interface of length L of the two squares consist of 2L constraints on the 8L 
parameters of the original two solution spaces. 

Remark 3.7 It is important to note that the neighbor subdomains may have different 
4L as dimensions of solution spaces. Similarly, the number 4L does not imply that each 
edge of a square subdomain bears exactly L parameters. In that case Observation ^. 6\ still 
holds: Whatever number of parameters born on the common interface by the Dirichlet, 
or Neumann, or D&N data (see Observation ^. 1\) will be eliminated for both subdomains 
after merging, and thus subtracted from the sum of dimensions of the two solution spaces. 

3.2 The algorithm 

It should be noted that our new algorithm is almost indentical to that of [1] in structure 
and complexity; they differ in solution space and in the merging strategy. The complex- 
ities for the old and new differ by a constant multiple, with the new more efficient by a 
considerable factor. 

The main advantage of the new is its simplicity in discretization and merging, two 
major difficulties still being reckoned with in existing approaches. The total wave based 
new algorithm never solve any scattering problem in constructing the total wave solution 
spaces (TWSS) for subdomains and in merging the solution spaces. The new only deals 
with the total waves in subdomains until the solution space is obtained for the entire 
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D, whereas the old methods solve a local scattering problem every step of the way, and 
have to deal with the typical corner singularities of the local scattered waves at the four 
corners of a subdomain. As two scattered wave solution spaces merge, four of the eight 
singularities at the eight corners of the two subdomains cancel out and disappear. This 
cancellation would give rise to conditioning problems. The total waves in a subdomain 
have natural extensions to the outside. They are unaware of the manmade corners of 
the subdomains, and are easier to sample and merge. 

We reiterate that for simplicity in describing our new algorithm, we assume that the 
scatterer q is a smooth function in D, and vanishes smoothly outside square D. In actual 
implementation, we will relax this requirement to include piecewise smooth scatterers 
with jumps inside and on the boundary of a fairly arbitrary, bounded domain D in 
two dimensions. We will also deal with jump discontinuities which form corners in the 
original medium; see §3.31 for more details. 

The total wave based fast direct solver is described in the following four subsections. 



3.2.1 Construct TWSS for bottom level subdomains 

Step 1 of the algorithm is to construct total wave solution spaces TWSS for each square 
subdomains S at the bottom level of the quadtree. 

TWSS is obtained by solving the homogeneous, variable coefficient Helmholtz equa- 
tion fl3T|) in S, subject to no boundary conditions. A high order method is provided 
in §3.3[ and also implemented numerically. Here we present a second order method to 
illustrate the procedure. The square S is discretized S a m-by-m uniform mesh Q. Five 
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Figure 1: Uniform mesh on a square S 



point stencil, to replace the Laplacian of (131) on every point, including those on the 
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boundary dS, produces m 2 homogeneous equations 



Ahu(x) + k 2 n(x)u(x) = 0, x E Q 



(4) 



The equations for the 4(m — 1) boundary points of the mesh will require 4m additional, 
free variables u(x), with x a step h away from the boundary; see Figured] for points x 
outside S marked by O. All together, there are m 2 equations for m 2 + 4m unknowns 
u(x), and so the null space of the discrete operator + k 2 n(x) is 4m dimensional, 
namely there are 4m nontrivial solutions to the m 2 equations. 

These 4m basis functions for the TWSS in S are collected in a matrix U s of size 
m 2 -by-4m. The basis functions are also evaluate at the boundary dS, and when paired 
with their normal derivatives, provide the D&N data for the TWSS. For simplicity, we 
assume that dS is sampled with 4m points, and that there are m points on each of the 
four edges of S, so that the matrices 



U = U„ 



as 



u 



(4m) 



], U n = {d n U s }\ ds = [u% 



(1) „(2) 



U 



(4m) 



are each of dimension 4m-by-4m. Let 



G 



U 
U„ 



G 



-U 



Un 



(5) 



(6) 



so G is of size 8m-by-4m. Our numerical experiment shows that these 4m bound- 
ary points can be equispaced, as opposed to crowded toward the four corners. Denser 
sampling points are required near a singularity of the solution (arising from medium 
discontinuity), or near locations where the total waves have more evanescent modes due 
to medium complexity. 



3.2.2 Merging two subdomains 

Step 2 of the algorithm is bottom up merging: On each level of the quadtree, merge two 
neighbor square subdomains, whose TWSS's are available, to construct TWSS for the 
resulting rectangular domain. Merge again two neighbor rectangles to form TWSS for 
the resulting square domain on the higher level. Stop at the highest level which contains 
only one square domain that is D. Merging is achieved by imposing continuity conditions 
on the D&N data on the common interface of two subdomains; see Observation 13.61 Let 
G\ and G2 be the D&N data matrix of the two square subdomains Si and S2. Let T be 
their common interface; see Figure [2j Then merging G\ and G2 to produce D&N data 
matrix G for S = S\ U S2 requires enforcing the continuity of the D&N data across T, 
by solution of the homogeneous linear system of 2m equations 

Gi| r n = -G 2 \r r 2 , namely [Gi\ r , G 2 \r] 



n 

r-2 



(7) 
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Figure 2: Merge two squares S\ and £2 



where Gi\r is the D&N data matrix Gi restricted on T; it is a matrix of size 2m-by-4m. 
The coefficient matrix [Gi|r, 62 1 r] is 2m-by-8m with a null space of dimension 6m. Let 
T\ and T 2 be matrices of size 4m-by-6m whose columns consist of the 6m solutions T\ 
and r 2 of ©. The D&N data matrix G on = {dSi \ T} U {<9S 2 \ T} are given by 

^l^sAr} = Gi|{ 9Sl \r} T u and G|{as 2 \r} = G 2 \{ds 2 \v} T 2 (8) 

We refer to (J7J), (JSJ) as the merging formulas, to be used throughout the bottom-up 
merging process. 

In the remainder of this subsection we define and determine splitting. Let u be a 
total wave - solution of the Helmholtz equation fl3Tj) in S = Si U S 2 . Let g, g 1 , g 2 be the 
D&N data for u on dS, dSi, dS 2 , respectively. Thus, there exist coefficients 7, 71, 72 
such that 

9 = G>y, g x = G171, g 2 = G 2 >y 2 (9) 



Definition 3.8 Splitting is the operation to determine gi,g 2 from g, in terms of their 



coefficients. The linear map S p : 7 1— > 



Ti 

72 



is referred to as the splitting operator. 



It follows OH]) immediately that S p , of size 8m-by-6m, is given by the formula 



S n 



T-2 



(10) 



Remark 3.9 All the matrix operations after Step 1 are carried out on the boundaries 
of the subdomains each with some 4m points on the boundary, instead of inside the 
subdomains each with about m 2 points. Therefore, the matrix operations are not nearly 
as costly. The entire merging step will cost only 0(k 3 ) flops for a k-by-k wavelength 
problem on D; see JJ^ for a complete analysis. 
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Subsections 13.2.11 and 13 . 2 . 21 are about constructing null spaces, and no scattering problem 
has been solved so far. 



3.2.3 Decompose TWSS via Green's formula 



Step 3 of the algorithm is to decompose TWSS for the whole scatterer D into the incoming 
and outgoing parts. Now that the TWSS is constructed for D with the D&N data matrix 
G available on dD which lies in the free space, the Green's identities apply, with the free 
space Green's function. In particular, the D&N data matrix G for the total waves can 
be split by the projector P_, see §4.21 for definition and technical details, to obtain the 
D&N data matrix Go for the incident parts of the total waves. 
Let the incident components of G be denoted by 

d n U 

Expressing the D&N data of the prescribed incident wave Uq on dD by the basis Go 



Gq 



so that G = P- G 



(11) 



G 7 



u 
d n u 



we solve (fT2l) to obtain the coefficients 7. Obviously, 



u 
d n u 



G 7 



(12) 



(13) 



is the D&N data of the total wave corresponding to the prescribed incident wave uq. For 
many applications, we want the scattered wave, and perhaps also its normal derivative, 
on dD. In this case, the D&N data of the scattered wave is obtained by subtracting 
(TI2"]) from ({TBI) . This concludes our algorithm if the scattered wave on dD is all we want; 
otherwise continue to the next step. 



3.2.4 Split a total wave 

Step 4 of the algorithm construct the total wave inside D. It is a top-down splitting 
process to propagate recursively along the quadtree the coefficient 7 from a domain S 
to its sub domains Si and £2. 

This is the reverse of Step 2 detailed in §3.2.2} see (TTOT) for details. Continue the 
recursive splitting till 7 is available for every bottom level square S. According to (JSJ), 
the total wave u in S is obtained by 

u = U sl (14) 

Now the total wave u of (1131) is available everywhere inside D. 

This is the end of our algorithm, and the volume scattering problem for a given 
incident wave uq is solved. 
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3.3 Some implementation details and remarks 

There are several very accurate methods to construct TWSS for a small subdomain S 
on the bottom level of the quadtree. We will present two typical methods: Collocation 
and the weak formulation. 



3.3.1 Collocation method to construct TWSS 

The solutions u of ( l3Tj) in S can be approximated by polynomials or bandlimited functions 
or some other suitable basis Bj(x). Thus we represent the total waves u in S by 

u(x) = CjBj(x), xES (15) 
j 

The homogeneous, variable coefficient Helmholtz equation fl3TT) . with u given above, is 
evaluated at some n suitable locations in S, giving rise to n homogeneous linear equations 
for cj 

Ac = (16) 

The vectors c from the null space of matrix A are then used to construct the total waves 
u for the TWSS in S. 

In numerical experiments, we used bandlimited functions of the form 

B(k, 8, x) = exp(ik(x 1 cos(#) + x 2 sm(9))), x e S, 6 e [0,2vr), ke[k 1 ,k 2 \ (17) 

The n collocation points on S are either uniform mesh - see Figured] - or a graded mesh 
such as the tensor legendre points. 



3.3.2 Weak formulation for TWSS 

The variational formulation for the Helmoltz equation (|3T|) offers more flexibility for an 
irregular subdomain S, and leads to typical finite element solution for ( !3T|) . Without 
boundary conditions for ( I3T1) . the weak formulation assumes the form 

-/v.V* + /^ + / tWlWl )^0 (18) 

JS JdS Js 

where u and v are both from a function space such as the one spanned by basis in (fT5]) . 
or by a typical finite element basis. In numerical experiments, we used ffTTj) as basis. 
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3.3.3 Sampling and completeness of TWSS 

For the prescribed precision e, the completeness of TWSS depends only on sampling, 
specifically the number of points and their distribution in each bottom level subdomain. 
Sampling rate is determined, as is well known, by the local wave number, by the complex- 
ity of the local medium, and by the distance to the nearest singularities arising from (i) 
the sources of incident wave (ii) the corners (iii) the edges. These sources have different 
strength of singularities, but if e is small it is hardly necessary to treat them differently 
in sampling rate near them. In the standard case when the medium changes smoothly 
and slowly relative to wavelength, and if the subdomain is far from a singularity, 2 to 4 
points per local wavelength is usually sufficient. 

On the other hand, since we always over sample, the TWSS will be unavoidably 
over complete to the prescribed precision. Once the TWSS space is constructed for 
a subdomain, particularly the bottom level ones, it may be compressed by SVD, or by 
pivoted Gram-Schmidt or QR. This will only make the TWSS healthier (because the total 
waves will be orthonormalized) . It will not reduce the level of over-completeness (because 
smaller singular values of matrix G of ^ may not correspond to more evanescent waves). 



3.3.4 Discontinuities in the medium 

If a bottom level subdomain S contains a smooth interface across which the index of 
refraction jumps, then S must be divided into two subdomains. These subdomans of S 
becoms bottom level subdomains. 

Let's assume now that a bottom level subdomain S contains a corner of the medium 
over which the index of refraction jumps. Then S must be partitioned into two subdo- 
mains, one contains the corner, the other is the complement. TWSS for each subdomain 
will contain some regular solutions and some singular solutions because of the corner. 
The regular ones will be constructed, say, by the collocation method ()15p . The singular 
ones will also be constructed by f lT5|) . except that the basis functions for them must be 
chosen to contain the local singular behavior of the total wave near the corner. Finally, 
the regular and singular solutions are collected together on the boundary and compressed 
by SVD or QR. 

As is well established, near such a corner the solution is spanned by the Bessel-Fourier 
terms of fractional orders, 



where x near the corner is assigned a polar coordinates (r, 6) centered at the corner. 




(19) 



V 



14 



3.3.5 Equations v.s. solutions 

The existing numerical methods for scattering or general elliptic problems can be divided 
into three categories, according to how much they are involved in building solution space. 

1. Our total wave approach deals with the total waves. The differential equations are 
treated only at the bottom level subdomains, and without boundary conditions. 
Merging TWSS of subdomains requires no knowledge of the PDEs. 

2. Methods based on scattering matrix or related objects [l]-[5] deal with the scattered 
waves for each subdomain on every level of the hierarchical domain decomposition 
(such as a quadtree). These methods solve a scattering problem on each bottom 
level subdomain, and they also need to know the underlying Green's function in 
merging two subdomains to construct the scattering matrix for the union of the 
two subdomains. Each subdomain on the hierarchy has manmade corners, which is 
visible to the scattered waves and show up as singularities in the scattering matrix. 

3. Finite difference, finite elements, or similar methods deals with the differential 
equations or their variational forms. Continuity conditions over an interface of 
subdomains are enforced on the equations rather than directly on the solutions, as 
our method does. 

Our total wave approach merges the subdomains directly and cleanly. Existing methods 
have to avoid dealing directly with the unwieldy analytical issues around the manmade 
corners arising from the manmade subdomains, by merging two "discretized" subdo- 
mains, or two overlapping ones. 

3.3.6 All merging formulas 

1. Merging two squares Si, S2 to a rectangle S. Let G\ and G2 be the two 

D&N data matrices to be merged to produce the D&N data matrix G for S = S\ U S2. 
Denote by Gij the part of Gj restricted to the edge shared with Gi, i ^ j. Then the 
merging-splitting matrices T\, T2 (see ((7j)) are solutions of the equation 

[G21 Gi2]2x8 

2. Merging two rectangles Si, S 2 to a square S. Let G\ and G 2 be the two 

D&N data matrices to be merged to produce the D&N data matrix G for S = Si U S2. 
Denote by G^ the part of Gj restricted to the edge shared with Gj, % ^ j. Then the 
merging-splitting matrices T\, T 2 are solutions of the equation 

[G2I Gi2]4xl2 



Tl 
T 2 



>2x6 



Sxfi 



(20) 



T-i 
T 2 



12x8 



(21) 
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3. Merging four squares Si, i = 1 : 4 to a square S. Let Gi, i = 1 : 4 be the four 
D&N data matrices to be merged to produce the D&N data matrix G for S = Uf =1 Si. 
Denote by Gij the part of Gj restricted to the edge shared with Gi, i j£ j. Then the 
merging-splitting matrices 7$, 1 — 1:4 are solutions of the equation 



G21 


G12 










" T x ' 





G32 


G23 







T 2 








G43 


G34 




T 3 


G41 








G\4 _ 


8x16 


T4 



3.4 Extensions to surface scattering problems 

Our total wave approach for VSP can be extended to a surface scattering problem (SSP). 
For simplicity, we consider a 2-D SSP off a sound soft (zero Dirichlet for total wave) 
smooth convex body D such as the unit disc with the sources of the incident wave uq 
well separated from D. Given uq on 3D, the SSP is to determine d n v, the normal 
derivative of the scattered wave v on 3D. 

The following steps outline a possible extension of the TWSS method to SSP. 

1. The bottom level subdomains. Dividing the annulus 1 < r < 1 + h, for some h > 0, 
along the radial direction into sufficient many pieces. Each piece 

Ai = {(r, 9), l<r<l + h, 0,_i < 9 < 9 t }, i = 1 : n (23) 

is bounded by four curves: two arcs and two straight radial segments. Remove the 
arc of r = 1 + h. The remaining three curves form the 2-th sub domain Tj on the 
bottom level: 

r fl = {(1,9), 9 1 ^ 1 <9<9 1 } (24) 
T l2 = {(r,0i_i), 1 <r < 1 + h} (25) 
T i3 = {(r,0 t ), l<r<l + h} (26) 

Ti = r a ur i2 ur i3 (27) 



2. TWSS in bottom level subdomains. For a prescribed precision construct TWSS in 
Ai subject to the zero Dirichlet condition (sound soft) on the arc Fa. Let there be 
m total waves in TWSS of the form = Uq ^ + j = 1 : m. For each the 
TWSS will contain four functions (i) and d n u^ restricted on the two radial 
segments rj 2 ,Fj 3 (ii) Uq and d n v^ on the arc r 4 i. 

3. Merging. Merge the subdomains, starting from the bottom level ones Fi, recur- 
sively and upward along a hierarchical domain decomposition. Two neighboring 



16 



subdomains C±, C2, separated by a common radial segment, will be merged to form 
their parent submain P by continuities of total wave and its normal derivative on 
the interface. 



4. After merging. The TWSS will contain pairs of functions on dD, the unit circle. 
Each pair is of the form 

uf,d n v®, 1=1: N (28) 

where d n v^' is the normal derivative of the scattered wave off D induced by an 
incident wave Uq , N is proportional to arclength of dD measured in wavelength. 

5. Construct d n v. Finally we construct the normal derivative, on dD, of the scattered 
wave v off D induced by the prescribed incident wave u$. Spanning uq by u 

u = J2^o ] (29) 

t 

we use the coefficients to produce d n v 

d n v = J2^d n v^ (30) 



This is the end of the algorithm. We have constructed the pair (v, d n v) on dD, with 
v = —Uq, and solved the SSP for sound soft scatterer D. The solution of (|29|) for ae is a 
bottleneck of the procedure. 



4 The Volume Scattering Problem 

One of the misfortunes of the 20th Century applied mathematics is that the volume 
scattering problem (VSP) for the Helmholtz equation was posed and is still being solved 
today as a boundary value problem for the scattered wave in a domain D, and worse yet, 
for each subdomain of D in a domain decomposition setting. The scattering problem is 
a very special boundary value problem in that both the Dirichlet and Neumann data of 
the incident wave are available on the boundary of the whole scatterer D. In this section 
we present the classical formulations for VSP, and special properties useful for the total 
wave based fast direct solver. 

Given index of refraction n(x) = yT-j- q(x) in a bounded domain D, we consider the 
volume scattering problem in /c-space governed by the Helmholtz equation 

Au + k 2 n 2 (x)u = 0, or An + k 2 (l + q)u = (31) 
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where u is the total wave, q is the scatterer, n = 1 and q = outside D. In a typical 
setting, m is of the form 

u = uq + (32) 

where the incident wave Uq is given in D and the scattered wave v : K 2 y C is to be 
determined as a solution of the inhomogeneous Helmholtz equation, 

Av + k 2 v = -k 2 q{u + v) (33) 

subject to the Sommerfeld radiation condition 

( dv \ 

lim yft — - ikv ) = (34) 

The scattering problem fl33l) . fl34|) can also be formulated as the Lippmann-Schwinger 
equation 

a(x) + k 2 q(x) [ G(x, = -k 2 q(x)u (x) (35) 

for the monopole density a in D, which is related to v by 

v(x)= [ G(x,t)a(t)dt (36) 
where G = —(i/4:)H (k\x — £|) is the free space Green's function. 



4.1 Dirichlet and Neumann Data for uo on dD 

For the scattering problem ( |33l and (J34l) . or ( |35l) and (J36l) . the incident wave uo must 
be available inside the scatterer D, or better, its sources outside D are prescribed. 

Observation 4.1 For a bounded domain D with a regular boundary 3D, the Dirichlet 
and Neumann data for the incident wave uq are always available on dD. 

To verify this statement, we observe that more often then not in a typical application, 
the incident wave u of a volume scattering problem is specified by its sources outside 
D, such as a monopole or a plane wave. In that case, the D&N data Uq and d n Uo can be 
evaluated directly on the boundary. 

Suppose Mo is only given in D, as required by ( 1331) or ( 1351) . By Green's third identity, 
Mo, being solution of Av + k 2 v = in D, is given by its D&N data 

u (x) = ~ J QD (dnMC) ■ G(x, - «o(0 • ds &, (37) 

So the D&N data for u can be recovered by solving f l37|) as an equation. This solution 
process can often be simplified if u is given not only in D, but also on D. 



18 



4.2 Green's Third Identity as the Interior Projector P 

Let W = L 2 (dD) x L 2 (dD). Let W± be two subspaces of W defined by the formulae 



W_ = { 0, d n v) | Av + k 2 v = in D } 

W+ = { (v, d n v) | Av + k 2 v = outside D subject to fl3l) } 



(38) 
(39) 



Therefore, W_ consists of the D&N data for incident waves in D, and W+ consists of the 
D&N data for outgoing (or scattered) waves outside D. 

For a bounded domain D with a regular boundary dD, the linear map (0, ip) £ i— > 
v £ L 2 (D) defined by Green's third identity 



v(x) 



dD 



dn{£) 



(40) 



produces v with Av + k 2 v = in D. Let x approach dD from .D, and denote the limit 
by x_ corresponding to x £ We thus obtain the linear map P_ : (0, ^) £ W i— )■ 



<9„w(x) 



so 



913 



a„0(O-G(x„,o-0(O 

<9G(z_,£) 



a„0(O 



dn(x) 



d 2 G(x.,Q 
dn{x)dn{£) 



ds(C), 



(41) 
(42) 



Likewise, using Green's third identity outside D, we introduce another linear map P + : 
(<f>,ijj) £ W H- (f,(9 n f) £ W+. In terms of the standard layer potential operators 
S, K, K' , T defined by (1351) - (1391) . the two operators P± are given by 



" " 




'10' 




' -K S ' 




" " 




-{i 


/ 


+ 


— T if' 


} 




" " 




"JO" 




" -Ji S 




" " 






I 




-T K' 


} 





(43) 
(44) 



It follows immediately from Green's third identity (1371) that 



Observation 4.2 TTie linear map P± :W ^ W± is a projector converting an arbitrary 
pair of boundary data (0, ip) £ W to the D&N data for an outgoing/incident wave. In 
particular, P_ maps the D&N data of the total wave u of UJty) to those of its incident 
component u , and P + maps the D&N data of total wave u to those of its scattered 
component v 
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Remark 4.3 Step 3 of the algorithm splits the D&N data of a total wave on dD for its 
scattered and incident parts using P_. Since the D&N data is nearly twice redundant, we 
only need the D&N data for the incident part evaluated at about half as many discretiza- 
tion points (quadrature nodes for P-). We will select these points away from corner and 
edges to avoid singularities of P_ there. 

In addition, if the wave number k is not a Dirichlet eigenvalue of the Laplacian in D, 
namely if k does not hit a geometric resonance of D, then only the first half of Go in 
(fT2l) is required to determine 7. In other words, the second half of P_, which involves 
the hyper singular kernel T, is not necessary for splitting. 



5 Layer potential representation for P + , P_ 

For x close to the the boundary dD, the single and double layer potentials 



p(x) 
q(x) 



G(x, 0^(0^(0, 



dD 



dG(x,£) 



1 3d 9n(0 

associated with arbitrary pair of functions (<p,ip), can be rewritten 



(45) 
(46) 



p(x ± hn(x)) 
q(x ± hn(x)) 



G(x±hn(x),0i>(0ds(t), 



dD 



dD 



8G(x ± hn{x),£) 



0(0^(0 



(47) 
(48) 



where x now is on the boundary dD. Taking the directional derivative of p, q in the 
normal direction n(x), we have 



dp(x ± hn(x)) 

dn(x) 
dq(x ± hn(x)) 



dG{x±hn(x),£) 

dn(x) 
d 2 G{x±hn{x),i) 



mds(0, 
mds(0 



dn(x) J dD dn(x)dn(^) 

Therefore, the Dirichlet and Neumann data (0, n ) of ( 140]) are given by 



p(x ) — q(x ) 
(p(x~) — q(x~)) 



(49) 
(50) 



' " 


= P 


' " 
.4. 


= — lim 


. 0™ . 







dn(x) 



(51) 
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for x = x — hn(x), x G dD. The Dirichlet and Neumann data 
by 

p(x + ) — q(x + 



i> n ) of (142]) are given 



' " 




" " 


= lim 






/i->+0 



3^y(p(x+)- 9 (x + )) 
for x + = x + hn(x), x G <9.D. The use of the jump conditions 



lim g(x ± hnix)) 



lim 



reduces (I5"T|) and 



<9p(x ± hn(x)) 
dn(x) 

to expressions for P± 



J 
/ 



± 



-K S 
-T K' 



where the layer potential operaters S, K, K', T are defined by 



(5a) (x) 
(Ka)(x) 
(K'a)(x) 

(Ta)(x) 



dD 

dG(x,Q 

dG(x,Q 
dn(x) 

d 2 G(x±hn(x),0 

hm 

h->+0 



(52) 

(53) 
(54) 

(55) 

(56) 
(57) 
(58) 
(59) 



f dD dn{x)dn{£) 

for x G 3D. Obviously, P+ + P- = I implying that the decomposition, of an arbitrary 
pair of data (<f), i/j) on boundary into incoming and outgoing parts, is complete. 

Remark 5.1 It is well-known that for a smooth dD the operators S, K, K' are bounded 
from C{dD) to C(dD), whereas T are bounded from C 1 (dD) to C(dD). 



Define the four linear operators S±, K±, K' ± , T± by the formulae 

(S ± ^)(x) 
(K ±( f>)(x) 



lim p(x ± hn(x)), 



(K' ± ^)(x) 



lim q(x ± hn(x)), 

dpix ± hnix)) 
hm — — 

h^+o onyx) 

.. dqix ± hnix)) 
hm 

h->+o dn[x) 



(60) 
(61) 

(62) 
(63) 
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for x G dD, and we see that 



K ± = K T \l, K' ± = K'±h, S ± = S, T ± = T. (64) 

As is well-known, the exterior Green's formula maps to zero the Dirichlet and Neumann 
data (0, <j) n ) of a outgoing wave <fi in D. Conversely, the interior Green's formula maps to 
zero the Dirichlet and Neumann data (ip, ip n ) of a outgoing wave ip outside D; therefore, 

P + oP^=P_oP + = 0, W„ n W+ = {0}. (65) 



Lemma 5.2 P is a projector if and only if there exists a unique operator Q such that 

1 , - 1 



Q 2 = -I, P=-I + Q 



(66) 



In particular, 



P± = -I± Q, with Q 



-K S 
-T K' 



(67) 



Furthermore, Q 2 = 1/4 implies that 



ST = --I + K 2 , TS = —-I + (K') 2 , KS = SK', TK = K'T (68) 

Finally, the Dirichlet-to- Neumann maps A± : <fi i— > <p n , with (0, n ) G W±(dD), are given 
by the formulae 



A- 



(69) 
(70) 



References 

[1] Y. Chen A fast direct solver for the Lippmann-Schwinger equation, Advances in 
Computational Mathematics, 16: 175-190, 2002. 

[2] P.G. Martinsson and V. Rokhlin A fast direct solver for boundary integral equations 
in two dimensions, Journal of computational physics, 205 (2005) Pages 1-23 



22 



[3] Y. Chen Rapid perturbational calculations for the Helmholtz equation in two dimen- 
sions. Discrete and Continuous Dynamical Systems - Series A (DCDS-A), Volume: 
18, Number: 4, August 2007, page 627 - 636 

[4] L. Gurel, W. Chew Fast direct (noniterative) solvers for integral- equation formula- 
tions of scattering problems, Antennas: Gateways to the Global Network, Vol. 1, 
IEEE Antennas and Propagation Society International Symposium, Vol. 1, IEEE 
Press, New York, 1998, pp. 298-301 

[5] W. Hackbusch A sparse matrix arithmetic based on H-matrices. Part I: Introduction 
to H-matrices, Computing 62. (1999), 89-108 

[6] E. Michielssen, A. Boag A multilevel matrix decomposition algorithm for analysing 
scattering from large structures, IEEE Trans. Antennas and Propagation 44 (8) 
(1996) 1086-1093. 

[7] L. Grasedyck, R. Kriemann, S. Le Borne Domain-decomposition based H-matrix 
preconditioners, Proceedings of DD16. LNSCE, vol. 55, 661-668. Springer, Berlin 
(2006) 

[8] S. Chandrasekaran, M. Gu, T. Pals A Fast ULV Decomposition Solver for Hierarchi- 
cally Semiseparable Representations, SIAM Journal on Matrix Anal. Appl, Volume 
28, Issue 3 (2006), 603 - 622 

[9] Y. Chen Riccati equations for scattering matrices on level surfaces, Inverse Prob- 
lems, vol 21, Pages 1745-1756, 2005 

[10] P.G. Martinsson, V. Rokhlin, and M. Tygert, On interpolation and integration in 
finite- dimensional spaces of bounded functions, Communications in Applied Math- 
ematics and Computational Science, 1, Jan. 2006. 



23 



